############################################################################
########          Keels, Joo, Wiegand  Replication Codes            ########
########          Models were run on a Windows system PC            ########
########  R version 4.2.1 (2022-06-23 ucrt) -- "Funny-Looking Kid"  ########
############################################################################

## Install Packages ## 

install.packages(foreign,logr, reshape2, forcats, dplyr, ggplot2)


#######################
##### Main Paper ######
#######################

library(foreign)
library(logr)
library(reshape2)
library(forcats)
library(dplyr)
library(ggplot2)


# set working directory 
wd <- setwd("C:/Users/MinHyung_Joo/Dropbox/ICW Peace Agreement paper/AJPS R&R/KJW_Replication Data")
wd <- getwd()

# read in data 
base <- read.dta("KJW_peaceduration_maindata.dta")

#################
## Figures 1-2 ##
#################

par(mar=c(4,4,1,1))

# plot and save as pdf files

pdf("Fig1a.pdf", width=6,height=6, paper='special')
hist(base$total_unad_pct, main="", cex=1.2, xlab="Proportion of Unaddressed Issues")
dev.off()

pdf("Fig2a.pdf", width=6,height=6, paper='special')
hist(base$intangible_unad_pct, main="", cex=1.2, xlab="Proportion of Unaddressed Intangible Issues", breaks = 10)
dev.off()

pdf("Fig1b.pdf", width=6,height=6, paper='special')
hist(base$total_unad_durmaj_pct, main="", cex=1.2, xlab="Proportion of Unaddressed Enduring Issues", breaks = 10)
dev.off()

pdf("Fig2b.pdf", width=6,height=6, paper='special')
hist(base$intang2_unad_durmaj_pct, main="", cex=1.2, xlab="Proportion of Unaddressed Enduring Intangible Issues", breaks = 10)
dev.off()

##############
## Figure 3 ## 
##############

# define analysis time frame
time <- seq(from=min(base$ltime), to=max(base$ltime), by=0.1)

### 20 % ###
# create blank matrix
nph <- matrix(NA, ncol=4, nrow=length(time))

# esimate hazard ratios
for(i in 1:length(time)){
  nph[i,1] <- exp(3.080/5) # Unaddressed 
  nph[i,2] <- exp((7.666-3.401*time[i])/5) # Enduring Unaddressed
  nph[i,3] <- exp(2.223/5)   # Intangible Unaddressed
  nph[i,4] <- exp((5.522-2.428*time[i])/5)   # Enduring Intangible Unaddressed
}

# name variables and reshape data
nph <- as.data.frame(nph)
names(nph) <- c("Unaddressed (long dash)", "Enduring Unaddressed", "Intangible Unaddressed", 
                "Enduring Intangible Unaddressed")

dd = melt(nph)
dd$time <- rep(time, 4)
names(dd) <- c("Key", "value", "time")

# plot and save as a pdf file
pdf("Fig3-20pct.pdf", width=6,height=6, paper='special')

ggplot(dd) +
  geom_line(aes(x=time, y=value, linetype=Key), size=0.8) +
  theme(panel.background = element_blank(), legend.position = "none",
        panel.border = element_rect(fill=NA, color="black")) +
  theme(axis.title.x  = element_text(size=12))+
  theme(axis.title.y  = element_text(size=12))+
  theme(axis.text.y=element_text(size=10),axis.text.x=element_text(size=10))+
  labs(y="Hazard Ratio", x="Log(time)") + 
  scale_linetype_manual(values=c("longdash", "solid", "twodash", "dotted")) + 
  geom_hline(yintercept=1, linetype='dashed', color='grey') + 
  annotate(geom="text", x=2.5, y=4.6, label="Unaddressed Rate: 20%")

dev.off()


### 40 % ###
# create blank matrix
nph <- matrix(NA, ncol=4, nrow=length(time))

# estimate hazard ratios 
for(i in 1:length(time)){
  nph[i,1] <- exp(3.080*.4) # Unaddressed 
  nph[i,2] <- exp((7.666-3.401*time[i])*.4) # Enduring Unaddressed
  nph[i,3] <- exp(2.223*.4)   # Intangible Unaddressed
  nph[i,4] <- exp((5.522-2.428*time[i])*.4)   # Enduring Intangible Unaddressed
}

# assign variable names and reshape data
nph <- as.data.frame(nph)
names(nph) <- c("Unaddressed (long dash)", "Enduring Unaddressed", "Intangible Unaddressed", 
                "Enduring Intangible Unaddressed")

dd = melt(nph)
dd$time <- rep(time, 4)
names(dd) <- c("Key", "value", "time")

# plot and save as pdf file 
pdf("Fig3-40pct.pdf", width=6,height=6, paper='special')

ggplot(dd) +
  geom_line(aes(x=time, y=value, linetype=Key), size=0.8) +
  theme(panel.background = element_blank(), legend.position ="none",
        panel.border = element_rect(fill=NA, color="black")) +
  theme(axis.title.x  = element_text(size=12))+
  theme(axis.title.y  = element_text(size=12))+
  theme(axis.text.y=element_text(size=10),axis.text.x=element_text(size=10))+
  labs(y="Hazard Ratio", x="Log(time)") + 
  scale_linetype_manual(values=c("longdash", "solid", "twodash", "dotted")) + 
  geom_hline(yintercept=1, linetype='dashed', color='grey')+
  annotate(geom="text", x=2.5, y=21, label="Unaddressed Rate: 40%")
dev.off()

### 60 % ###
# create blank matrix
nph <- matrix(NA, ncol=4, nrow=length(time))

# estimate hazard ratios 
for(i in 1:length(time)){
  nph[i,1] <- exp(3.080*.6) # Unaddressed 
  nph[i,2] <- exp((7.666-3.401*time[i])*.6) # Enduring Unaddressed
  nph[i,3] <- exp(2.223*.6)   # Intangible Unaddressed
  nph[i,4] <- exp((5.522-2.428*time[i])*.6)   # Enduring Intangible Unaddressed
}

# assign variable names and reshape data
nph <- as.data.frame(nph)
names(nph) <- c("Unaddressed (long dash)", "Enduring Unaddressed", "Intangible Unaddressed", 
                "Enduring Intangible Unaddressed")

dd = melt(nph)
dd$time <- rep(time, 4)
names(dd) <- c("Key", "value", "time")

# plot and save as pdf file 
pdf("Fig3-60pct.pdf", width=6,height=6, paper='special')

ggplot(dd) +
  geom_line(aes(x=time, y=value, linetype=Key), size=0.8) +
  theme(panel.background = element_blank(), legend.position = "none",
        panel.border = element_rect(fill=NA, color="black")) +
  theme(axis.title.x  = element_text(size=12))+
  theme(axis.title.y  = element_text(size=12))+
  theme(axis.text.y=element_text(size=10),axis.text.x=element_text(size=10))+
  labs(y="Hazard Ratio", x="Log(time)") + 
  scale_linetype_manual(values=c("longdash", "solid", "twodash", "dotted")) + 
  geom_hline(yintercept=1, linetype='dashed', color='grey')+
  annotate(geom="text", x=2.5, y=98, label="Unaddressed Rate: 60%")

dev.off()

### 80% ###
# create blank matrix
nph <- matrix(NA, ncol=4, nrow=length(time))

# estimate hazard ratios 
for(i in 1:length(time)){
  nph[i,1] <- exp(3.080*.8) # Unaddressed 
  nph[i,2] <- exp((7.666-3.401*time[i])*.8) # Enduring Unaddressed
  nph[i,3] <- exp(2.223*.8)   # Intangible Unaddressed
  nph[i,4] <- exp((5.522-2.428*time[i])*.8)   # Enduring Intangible Unaddressed
}

# assign variable names and reshape data
nph <- as.data.frame(nph)
names(nph) <- c("Unaddressed (long dash)", "Enduring Unaddressed", "Intangible Unaddressed", 
                "Enduring Intangible Unaddressed")

dd = melt(nph)
dd$time <- rep(time, 4)
names(dd) <- c("Key", "value", "time")

# plot and save as pdf file 
pdf("Fig3-80pct.pdf", width=6,height=6, paper='special')

ggplot(dd) +
  geom_line(aes(x=time, y=value, linetype=Key), size=0.8) +
  theme(panel.background = element_blank(), legend.position = "none",
        panel.border = element_rect(fill=NA, color="black")) +
  theme(axis.title.x  = element_text(size=12))+
  theme(axis.title.y  = element_text(size=12))+
  theme(axis.text.y=element_text(size=10),axis.text.x=element_text(size=10))+
  labs(y="Hazard Ratio", x="Log(time)") + 
  scale_linetype_manual(values=c("longdash", "solid", "twodash", "dotted")) + 
  geom_hline(yintercept=1, linetype='dashed', color='grey')+
  annotate(geom="text", x=2.5, y=450, label="Unaddressed Rate: 80%")

dev.off()


##############
## Figure 4 ##
##############

# assign variable names and structure data
varname <- c("Political", "Economic", "Natural", "Religious", "Ethnic", "Self-Determination", "Government")
coef <- c(9.213, -0.158, 1.375, 8.952, 1.277, 1.116, 0.941)
se <- c(1.934, 0.710, 0.991, 1.278, 0.710, 0.753, 0.816)
u95 <- coef + 1.96*se
l95 <- coef - 1.96*se

results <- data.frame(as.matrix(cbind(varname, coef, se, u95, l95)))
CIdata <- data.frame(model=varname, mean=coef, lower = l95, upper= u95)

# plot and save as pdf file 
pdf("Fig4.pdf", width=6,height=4, paper='special')

ggplot() +
  geom_errorbar(data=CIdata, mapping=aes(x=fct_rev(fct_inorder(model)), ymin=lower, ymax=upper), width=0.2, size=0.7, color="black")+
  theme_bw() + theme(panel.grid.major = element_blank(),
                     panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"),
                     axis.title.y = element_text(size=12, face="bold"),
                     axis.title.x = element_text(size=12, face="bold"),
                     axis.text.x = element_text(size=9, face='bold'),
                     axis.text.y = element_text(size=9, face='bold'))  +
  geom_point(data=CIdata, mapping=aes(x=model, y=mean), size=2,  fill="white") +
  geom_hline(yintercept=0, linetype="dashed",  size=0.5) + 
  xlab("Unaddressed Enduring Issue") + ylab("Coefficient") + 
  coord_flip()

dev.off()


##############
## Figure 5 ##
##############

#20%
# create blank matrix
nph2 <- matrix(NA, ncol=4, nrow=length(time))

# estimate hazard ratios 
for(i in 1:length(time)){
  nph2[i,1] <- exp((8.788 -3.978*time[i])*.2) # Unaddressed 
  nph2[i,2] <- exp((8.026-3.088*time[i])*.2) # Enduring Unaddressed
  nph2[i,3] <- exp((4.949-2.395*time[i])*.2)   # Intangible Unaddressed
  nph2[i,4] <- exp((3.973-1.410*time[i])*.2)   # Enduring Intangible Unaddressed
}

# assign variable names and reshape data
nph2 <- as.data.frame(nph2)
names(nph2) <- c("Unaddressed(long dash)", "Enduring Unaddressed", "Intangible Unaddressed", 
                 "Enduring Intangible Unaddressed")

dd2 = melt(nph2)
dd2$time <- rep(time, 4)
names(dd2) <- c("Key", "value", "time")


# plot and save as pdf file 
pdf("Fig5-20pct.pdf", width=6,height=6, paper='special')
ggplot(dd2) +
  geom_line(aes(x=time, y=value, linetype=Key), size=0.8) +
  theme(panel.background = element_blank(), legend.position = "none",
        panel.border = element_rect(fill=NA, color="black")) +
  theme(axis.title.x  = element_text(size=12))+
  theme(axis.title.y  = element_text(size=12))+
  theme(axis.text.y=element_text(size=10),axis.text.x=element_text(size=10))+
  labs(y="Hazard Ratio", x="Log(time)") + 
  scale_linetype_manual(values=c("longdash", "solid", "twodash", "dotted"))+ 
  geom_hline(yintercept=1, linetype='dashed', color='grey')+
  annotate(geom="text", x=2.5, y=6, label="Unaddressed Rate: 20%")

dev.off()

# 40%
# create blank matrix
nph2 <- matrix(NA, ncol=4, nrow=length(time))

# estimate hazard ratios 
for(i in 1:length(time)){
  nph2[i,1] <- exp((8.788 -3.978*time[i])*.4) # Unaddressed 
  nph2[i,2] <- exp((8.026-3.088*time[i])*.4) # Enduring Unaddressed
  nph2[i,3] <- exp((4.949-2.395*time[i])*.4)   # Intangible Unaddressed
  nph2[i,4] <- exp((3.973-1.410*time[i])*.4)   # Enduring Intangible Unaddressed
}

# assign variable names and reshape data
nph2 <- as.data.frame(nph2)
names(nph2) <- c("Unaddressed(long dash)", "Enduring Unaddressed", "Intangible Unaddressed", 
                 "Enduring Intangible Unaddressed")

dd2 = melt(nph2)
dd2$time <- rep(time, 4)
names(dd2) <- c("Key", "value", "time")


# plot and save as pdf file 
pdf("Fig5-40pct.pdf", width=6,height=6, paper='special')
ggplot(dd2) +
  geom_line(aes(x=time, y=value, linetype=Key), size=0.8) +
  theme(panel.background = element_blank(), legend.position = "none",
        panel.border = element_rect(fill=NA, color="black")) +
  theme(axis.title.x  = element_text(size=12))+
  theme(axis.title.y  = element_text(size=12))+
  theme(axis.text.y=element_text(size=10),axis.text.x=element_text(size=10))+
  labs(y="Hazard Ratio", x="Log(time)") + 
  scale_linetype_manual(values=c("longdash", "solid", "twodash", "dotted"))+ 
  geom_hline(yintercept=1, linetype='dashed', color='grey')+
  annotate(geom="text", x=2.5, y=33, label="Unaddressed Rate: 40%")
dev.off()

# 60%
# create blank matrix
nph2 <- matrix(NA, ncol=4, nrow=length(time))

# estimate hazard ratios 
for(i in 1:length(time)){
  nph2[i,1] <- exp((8.788 -3.978*time[i])*.6) # Unaddressed 
  nph2[i,2] <- exp((8.026-3.088*time[i])*.6) # Enduring Unaddressed
  nph2[i,3] <- exp((4.949-2.395*time[i])*.6)   # Intangible Unaddressed
  nph2[i,4] <- exp((3.973-1.410*time[i])*.6)   # Enduring Intangible Unaddressed
}

# assign variable names and reshape data
nph2 <- as.data.frame(nph2)
names(nph2) <- c("Unaddressed(long dash)", "Enduring Unaddressed", "Intangible Unaddressed", 
                 "Enduring Intangible Unaddressed")

dd2 = melt(nph2)
dd2$time <- rep(time, 4)
names(dd2) <- c("Key", "value", "time")

# plot and save as pdf file 
pdf("Fig5-60pct.pdf", width=6,height=6, paper='special')

ggplot(dd2) +
  geom_line(aes(x=time, y=value, linetype=Key), size=0.8) +
  theme(panel.background = element_blank(), legend.position = "none",
        panel.border = element_rect(fill=NA, color="black")) +
  theme(axis.title.x  = element_text(size=12))+
  theme(axis.title.y  = element_text(size=12))+
  theme(axis.text.y=element_text(size=10),axis.text.x=element_text(size=10))+
  labs(y="Hazard Ratio", x="Log(time)") + 
  scale_linetype_manual(values=c("longdash", "solid", "twodash", "dotted"))+ 
  geom_hline(yintercept=1, linetype='dashed', color='grey')+
  annotate(geom="text", x=2.5, y=195, label="Unaddressed Rate: 60%")
dev.off()


# 80%
# create blank matrix
nph2 <- matrix(NA, ncol=4, nrow=length(time))

# estimate hazard ratios 
for(i in 1:length(time)){
  nph2[i,1] <- exp((8.788 -3.978*time[i])*.8) # Unaddressed 
  nph2[i,2] <- exp((8.026-3.088*time[i])*.8) # Enduring Unaddressed
  nph2[i,3] <- exp((4.949-2.395*time[i])*.8)   # Intangible Unaddressed
  nph2[i,4] <- exp((3.973-1.410*time[i])*.8)   # Enduring Intangible Unaddressed
}

# assign variable names and reshape data
nph2 <- as.data.frame(nph2)
names(nph2) <- c("Unaddressed(long dash)", "Enduring Unaddressed", "Intangible Unaddressed", 
                 "Enduring Intangible Unaddressed")

dd2 = melt(nph2)
dd2$time <- rep(time, 4)
names(dd2) <- c("Key", "value", "time")

# plot and save as pdf file 
pdf("Fig5-80pct.pdf", width=6,height=6, paper='special')

ggplot(dd2) +
  geom_line(aes(x=time, y=value, linetype=Key), size=0.8) +
  theme(panel.background = element_blank(), legend.position = "none",
        panel.border = element_rect(fill=NA, color="black")) +
  theme(axis.title.x  = element_text(size=12))+
  theme(axis.title.y  = element_text(size=12))+
  theme(axis.text.y=element_text(size=10),axis.text.x=element_text(size=10))+
  labs(y="Hazard Ratio", x="Log(time)") + 
  scale_linetype_manual(values=c("longdash", "solid", "twodash", "dotted"))+ 
  geom_hline(yintercept=1, linetype='dashed', color='grey')+
  annotate(geom="text", x=2.5, y=1100, label="Unaddressed Rate: 80%")

dev.off()


###########################
##### Online Appendix #####
###########################

wd <- setwd("Your_Directory/KJW_Replication Data/")
wd <- getwd()

base <- read.dta("KJW_peaceduration_maindata.dta")


library(forcats)
library(dplyr)
library(ggplot2)

# assign variable names and structure data
varname <- c("Unaddressed-No Lootable", "Unaddressed-Lootable", "Enduring Unaddressed-No Lootable", "Enduring Unaddressed-Lootable")
coef <- c(2.464777, 4.722404, 2.642243, 5.377206)
se <- c(1.958459, 1.060295,  2.256953 , .9888235)
u95 <- coef + 1.96*se
l95 <- coef - 1.96*se

results <- data.frame(as.matrix(cbind(varname, coef, se, u95, l95)))
CIdata <- data.frame(model=varname, mean=coef, lower = l95, upper= u95)

# plot and save as pdf file 
pdf("FigA1.pdf", width=7,height=3.5, paper='special')

ggplot() +
  geom_errorbar(data=CIdata, mapping=aes(x=fct_rev(fct_inorder(model)), ymin=lower, ymax=upper), width=0.2, size=0.7, color="black")+
  theme_bw() + theme(panel.grid.major = element_blank(),
                     panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"),
                     axis.title.y = element_text(size=12, face="bold"),
                     axis.title.x = element_text(size=12, face="bold"),
                     axis.text.x = element_text(size=12, face='bold'),
                     axis.text.y = element_text(size=15, face='bold'))  +
  geom_point(data=CIdata, mapping=aes(x=model, y=mean), size=2,  fill="white") +
  geom_hline(yintercept=0, linetype="dashed",  size=0.5) + 
  xlab("") + ylab("Coefficient") + 
  coord_flip()

dev.off()

#### Intangible ####

# assign variable names and structure data
varname <- c("Intangible Unaddressed-No Lootable", "Intangible Unaddressed-Lootable", "Enduring Intangible Unadressed-No Lootable", "Enduring Intangible Unaddressed-Lootable")
coef <- c(   1.693538 , 2.252546,   1.502792 ,  2.894981)
se <- c( 1.555992  , .8424192,  1.726187 ,  .9756621)
u95 <- coef + 1.96*se
l95 <- coef - 1.96*se


results <- data.frame(as.matrix(cbind(varname, coef, se, u95, l95)))
CIdata <- data.frame(model=varname, mean=coef, lower = l95, upper= u95)

# plot and save as pdf file 
pdf("FigA2.pdf", width=6,height=3.5, paper='special')

ggplot() +
  geom_errorbar(data=CIdata, mapping=aes(x=fct_rev(fct_inorder(model)), ymin=lower, ymax=upper), width=0.2, size=0.7, color="black")+
  theme_bw() + theme(panel.grid.major = element_blank(),
                     panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"),
                     axis.title.y = element_text(size=12, face="bold"),
                     axis.title.x = element_text(size=12, face="bold"),
                     axis.text.x = element_text(size=12, face='bold'),
                     axis.text.y = element_text(size=15, face='bold'))  +
  geom_point(data=CIdata, mapping=aes(x=model, y=mean), size=2,  fill="white") +
  geom_hline(yintercept=0, linetype="dashed",  size=0.5) + 
  xlab("") + ylab("Coefficient") + 
  coord_flip()

dev.off()

### Ext Support ###

# assign variable names and structure data
varname <- c("Unaddressed-No External Support", "Unaddressed-External Support", "Enduring Unadressed-No External Support", "Enduring Unaddressed-External Support")
coef <- c( 41.23326,  -.7707167,  10.73296 ,  -6.972542)
se <- c(18.94516 , 6.439783 , 17.46793 ,  21.61242)
u95 <- coef + 1.96*se
l95 <- coef - 1.96*se


results <- data.frame(as.matrix(cbind(varname, coef, se, u95, l95)))
CIdata <- data.frame(model=varname, mean=coef, lower = l95, upper= u95)

# plot and save as pdf file 
pdf("FigA3.pdf", width=7,height=3.5, paper='special')

ggplot() +
  geom_errorbar(data=CIdata, mapping=aes(x=fct_rev(fct_inorder(model)), ymin=lower, ymax=upper), width=0.2, size=0.7, color="black")+
  theme_bw() + theme(panel.grid.major = element_blank(),
                     panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"),
                     axis.title.y = element_text(size=12, face="bold"),
                     axis.title.x = element_text(size=12, face="bold"),
                     axis.text.x = element_text(size=12, face='bold'),
                     axis.text.y = element_text(size=15, face='bold'))  +
  geom_point(data=CIdata, mapping=aes(x=model, y=mean), size=2,  fill="white") +
  geom_hline(yintercept=0, linetype="dashed",  size=0.5) + 
  xlab("") + ylab("Coefficient") + 
  coord_flip()

dev.off()



# End of file 

